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ABSTRACT 


A mathematical theory and an accompanying numerical scheme have been 
developed for predicting the oxidation behavior of carbon silicon carbide (C/SiC) 
composite structures. The theory is derived from the mechanics of the flow of ideal gases 
through a porous solid. The result of the theoretical formulation is a set of two coupled 
nonlinear differential equations written in terms of the oxidant and oxide partial 
pressures. The differential equations are solved simultaneously to obtain the partial vapor 
pressures of the oxidant and oxides as a function of the spatial location and time. The 
local rate of carbon oxidation is determined using the map of the local oxidant partial 
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vapor pressure along with the Arrhenius rate equation. The nonlinear differential 
equations are cast into matrix equations by applying the Bubnov-Galerkin weighted 
residual method, allowing for the solution of the differential equations numerically. The 
numerical method is demonstrated by utilizing the method to model the carbon oxidation 
and weight loss behavior of C/SiC specimens during thermogravimetric experiments. The 
numerical method is used to study the physics of carbon oxidation in carbon silicon 
carbide composites. 

Keywords: A. Carbon fibers; B. Oxidation; C. Modeling; D. Diffusion, Porosity; 

1. INTRODUCTION 

The ability of carbon fiber-reinforced silicon carbide composites (C/SiC) to 
maintain its strength and stiffness at high temperatures as well as its low density make it 
an attractive candidate for many applications in future spacecraft. These applications 
include turbomachinery components and thrust chambers in future propulsion systems as 
well as control surfaces, leading edges and thermal protection systems for vehicle 
airframes. One of the more formidable obstacles to the widespread use of C/SiC 
structures in future launch vehicles is that the carbon fibers oxidize at medium to high 
temperatures in an environment in which oxygen is present. This does not forbid the use 
of C/SiC in future launch vehicle applications, as long as it can be verified through 
testing and analysis that the component will maintain its strength and stiffness throughout 
its service life, with the demonstration of sufficient safety factors. As such, an assessment 
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of the oxidation behavior of C/SiC composite structures must be included along with the 
usual design analysis activities such as the thermal, dynamic and thermostructural 
analysis of the component. It is therefore necessary to develop a tool that is capable of 
determining the spatial distribution of the extent of oxidation and the residual strength 
and stiffness in the C/SiC component as a function of the time, temperature and 
environmental oxygen concentrations to which the C/SiC structure is exposed. Currently, 
no such oxidation analysis tool is available to designers, who wish to utilize C/SiC 
composites. 

Oxygen attacks the carbon in C/SiC composites both on the surface and in the 
interior of the composite. The oxygen achieves access to the interior of the composite via 
an interconnected network of passageways, which are formed by the combination of 
matrix cracks and void spaces, which exist both within the fiber bundles as well as 
between adjacent plies. Oxygen may also flow through separations of the fiber/coating 
and coating/matrix interfaces. The matrix cracks and interfacial separations are due to 
tensile stresses, which are a result of the thermal expansion mismatch between the carbon 
fibers and the silicon carbide matrix in concert with the temperature excursions during 
processing and cool down [1]. The large free spaces between plies and the void spaces 
within the fiber bundles are both due to insufficient void filling during matrix infiltration. 
The void spaces and matrix cracks are illustrated in Fig. 1 where the microstructure of a 
2-D C/SiC composite is shown. The oxidation of carbon in the interior of C/SiC 
composites is strictly tied to the transport of oxygen into and the transport of oxides out 
of the material. Any viable oxidation model for C/SiC composites must include the 
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solution of species transport equations as the transport has a direct impact on the rate of 
carbon oxidation. 

Oxidation models have been developed in the past in order to study the physics of 
the oxidation process in carbon fiber-reinforced composites. Medford [2] proposed a 
model to predict the oxidation of the Space Shuttle’s carbon-carbon wing leading edge, 
by simulating the diffusion of oxygen to the carbon-carbon substrate down a fissure in the 
SiC coating. Eckel, et al. [3] proposed a similar model to determine the oxidation 
recession rate of a single carbon fiber embedded in a non-reactive matrix. Halbig [4] 
adapted Eckel, et al.’s model to simulate the fiber surface recession in C/SiC specimens. 
His approach presupposes a crack extending through the specimen gage section, bridged 
by an array of continuous carbon fibers. All three models assume steady state diffusion of 
the oxidant to the site of oxidation. 

Although these previous studies have provided insight into the physics of carbon 
oxidation in ceramic composites, these approaches are not readily applicable to support 
the design analysis of C/SiC structures as they are impractical for predicting the residual 
strength and stiffness as a function of space and time for any arbitrarily-shaped C/SiC 
structure. Indeed, these previous methods study the problem of carbon oxidation on a 
very fine scale, retaining the heterogeneous nature of the composite. As a result, these 
methods would be too cumbersome if they were applied to analyze C/SiC structures on a 
global level. For this purpose, a numerical model that treats the C/SiC composite as a 
continuum is needed. 
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The purpose of this paper is to describe the development of an accurate analysis 
method that simulates the oxidation behavior of C/SiC composite structures in high 
temperature applications. In this paper, the mathematical foundation of the method is 
presented. The theory is derived by assuming that the C/SiC material is a homogeneous, 
orthotropic porous body with a solid skeleton that is a mixture of multiple solid 
constituents where some of the constituents are reactive. The oxygen (oxidant) and the 
oxides (product) flow through the pore network and the partial pressures of the gases vary 
with space and time. The pore volume of this porous body represents the collective 
volumes of the matrix cracks, the fiber/coating and coating/matrix interfacial separations, 
and both types of void volumes (inter-ply and intra-bundle). Applying the fundamentals 
of porous media to this problem, namely the mass conservation equation for each gas 
specie as well as the transport mechanisms, a set of coupled, nonlinear differential 
equations is obtained. The solution of these differential equations yields the partial vapor 
pressures of the oxidant and oxides as a function of space and time. The local rate of 
carbon oxidation is determined as a function of space and time using the map of the local 
oxidant partial vapor pressure along with the Arrhenius rate equation. The Bubnov- 
Galerkin weighted residual method is used to cast the governing differential equations 
into matrix equations in order to perform the solution of the differential equations 
numerically. 

In the application section of this paper, we demonstrate the use of the numerical 
method by applying it to simulate the carbon oxidation and weight loss behavior of C/SiC 
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specimens during thermo gravimetric analysis (TGA) experiments. Using the numerical 
method, we deduce the value of two key material parameters as well as the variation of 
these parameters with temperature and we demonstrate the profound influence of 
temperature on the oxidation behavior. 


2. THEORETICAL FORMULATION 

2.1 Mass Continuity for Flow Through Porous Media 

In order to develop the theory to model the oxidation process, we make use of the 
basic principals of porous media theory, namely the equation for the continuity of mass 
of gaseous species flowing through a porous solid as well as mechanisms for transport in 
the porous solid. It is assumed that the oxygen and the products of the oxidation reaction 
(oxides) exist only in the pores of the material in a gas form and that these gases behave 
ideally, that the solid skeleton consists of a mixture of both silicon carbide and carbon in 
the solid form and that the carbon is oxidized at the interface between the solid skeleton 
and the pore, in other words, at the wall of the pore space. 

The local form of the mass continuity equations for oxygen and the oxide species, 
flowing through a porous solid body, may be written as 

^ + < la > 
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and 


t^ + V J co=-h^K+*pc). (lb) 

respectively, where 0 is the volumetric porosity and V is the gradient operator and 

where pf is the local partial density in the pore volume, Mj is the molecular weight and 
is the local mass flux vector for specie i. The subscripts ox, co and c refer to the 
oxygen, the oxide (either carbon monoxide or carbon dioxide) and the carbon (solid 
form) species, respectively. 

In equations (la) and (lb), the symbols 91 and < R pc denote the local time rate 

of carbon fiber mass and pyrocarbon coating mass consumption due to the oxidation 
reaction per unit bulk volume, respectively. Further, and A 2 are the stoichiometric 
constants for the oxidation reactions. That is, is the ratio of the number of moles of 
oxygen consumed in the oxidation reaction to the number of moles of carbon consumed 
in the reaction and is the ratio of the number of moles of oxide produced in the 
oxidation reaction to the number of moles of carbon consumed. As such, the right hand 
side (RHS) of (la) is the local rate of oxygen mass consumed in the oxidation reaction 
per unit bulk volume and the RHS of (lb) is the local rate of oxide mass produced by the 
oxidation reaction per unit bulk volume. 
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Upon substituting the ideal gas law for each specie i (pf = PiM { / RT ) into the 
first term in equations (la) and (lb) and upon performing the differentiation, we obtain, 


after rearranging, 



(2a) 


and 



(2b) 


2.2 Mass Flux Constitutive Relations 

We will allow for the diffusion of oxygen and the oxides through the pore 
network via two mechanisms: gas pressure gradient-driven flow and concentration 
gradient-driven flow. Thus, we may write the mass flux vector of specie i as the sum of 
two mass flux vectors as 


J;=J 



( 3 ) 
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where Jf is the mass flux associated with the average velocity of the gas mixture and is 

gas pressure gradient-driven and where jf is the flux of specie / relative to the mixture 
average velocity and is concentration gradient-driven. 

The expression for the gas pressure gradient-driven flow is attributed to Darcy [5] 
and is written as 


— kVp (4) 

"g 

where k is the second-order material permeability tensor, ju g is the viscosity of the gas 
mixture and p is the total gas pressure. 

The concentration gradient-driven flow is given by modifying Fick’s law [6] as 

jf=-p p D AB < r'-vfe, (5) 

P P 

where p p is the local gas mixture density in the pore and D AB is the diffusivity of gas 

specie A with respect to gas specie B, and (p A is the second-order areal porosity tensor. 
The diffusivity D AB is given by the Chapman-Enskog equation which is derived from 
the kinetic theory of gases [6]. 
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The areal porosity tensor is a measure of the resistance to concentration gradient- 
driven flow through the pore network. We note that Bacos, et al. [7] use the ratio e I t to 
represent this resistance, where the quantity e denotes the porosity and t is defined as a 
tortuosity factor. In the present study, we will use the areal porosity tensor to represent 

the ratio eh. As such, <p A represents the combined resistance to concentration 
gradient-driven diffusion due to both the tortuosity and the fraction of the cross-sectional 
area which is occupied by pore volume. The value of the areal porosity will, of course, 
depend on the pore morphology. 

Making use of the ideal gas law for each specie i ( pf = PjMj / RT ) as well as the 
ideal gas law for the mixture ( p p = pM g / RT ), we can rewrite equation (5) in terms of 
the partial pressures p t : 


jf =-p P D AB -£ L <t A V( £L ), 

M g p 


where M g is the molecular weight of the gas mixture. 


Substituting equations (4) and (6) into (3) leads to 


( 6 ) 
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Jo, =-p&— k Vp OT -D AB ^jr—V A Vp ox 
Mg K1 P 

(7a) 

- p p' k .V Pm+DAB M^£*L <t A -VPco 
Mg K1 p 

and 

Jco =-pfo-l~k Vp co -D AB ^-^ A Vp co 
Mg K1 P 

(7b) 

Vp OT . 

// g Ai 

In the derivation of equations (7a) and (7b), it was necessary to employ Dalton’s law for 
the gas mixture (p = p ox + p co ), the distributive property of the gradient operator 

( Vp = Vp ox + Vp co ) and the ideal gas law for the gas mixture. 

It is apparent that, upon substituting equations (7a) and (7b) into (2a) and (2b), we 
will arrive at two coupled nonlinear differential equations written in terms of the partial 
pressures. These equations will be coupled, since terms involving the gradients Vp ox and 

VPco will appear in both equations. These equations will be nonlinear, since (7a) and 

(7b) contain terms that involve the product of the partial densities and the gradients of the 
partial pressures as well as terms that involve the product of the partial pressures and the 
gradients of the partial pressures. The two nonlinear differential equations must be solved 
simultaneously at each time step to obtain the partial pressures at each spatial location. In 


11 



the following sections, a numerical approach is presented to perform the simultaneous 
solution of these equations. 

2.3 Determination of Oxidation Reaction Rates 

The time rate of the carbon oxidation reaction is a function of the absolute 
temperature and the vapor pressure of the oxidant [8]. The dependence of the reaction 
rate on temperature and pressure is given by the Arrhenius rate equation. As we are 
concerned with a solid mixture containing carbon, then for any unit volume of material, 
we can write the Arrhenius equation in terms of a density ratio as 


where p c and p° are the instantaneous mass density and initial mass density of carbon 
in the solid mixture, respectively, and where k Q is the pre-exponential coefficient, E a is 
the activation energy of the oxidation reaction and n is the order of the reaction. The 
Arrhenius constants k Q , E a and n are obtained by curve fitting thermogravimetric 
analysis (TGA) measurements to the Arrhenius equation. 

In applying equation (8), it is necessary to recognize a few key differences 
between the oxidation of carbon fibers or pyrocarbon coating in TGA tests and the 
oxidation of these constituents in C/SiC composites. First, the ambient partial pressure of 


,( \ 
d Pj. 



( 8 ) 



P J 

\^c J oxidation 
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oxygen in the TGA experiments is constant and specified as a test condition, whereas the 
oxygen partial pressure in the C/SiC composite varies with time and spatial position. 
Second, the carbon surface area that is exposed to oxygen vapor may be approximated as 
close to 100% in the TGA experiments, whereas this surface area fraction in the C/SiC 
composite is considerably less than this estimate. It is therefore necessary to modify 
equation (8) in order to use the Arrhenius rate equation to determine the local rate of 
oxidation in C/SiC composites. 

We note that by definition, 91 c y = {^Pcf ^ t ) oxidation an d 

< R pc = {dp p C ! dt) oxidation , where p c f and p pc are the local mass density of carbon 

fiber and local mass density of pyrocarbon coating, respectively. The local mass densities 
and the local volume fractions are related by p c f = p c v c f and p pc = p c v pc , where v c y 

and v pc are the volume fraction of carbon fiber and pyrocarbon coating, respectively and 

where p c is the intrinsic density of carbon. 

Given these considerations, equation (8) may be rewritten separately for both the 
carbon fibers and pyrocarbon coating, as, 


% - ~Pc v °cf 


f \ 
Pox 

\P*0XJ 


n cf 


y/ c f kf expl 


'-if 

RT 


(9a) 


and 
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v = - 


Pc v pc 


( \ n pc 
Pox 
* 

\Pox ) 


y/ pc k pc exp 


'-j?' 

RT 


(9b) 


where and v pc are the initial values of the volume fractions and where k c J , E c f 

and n c j- are the Arrhenius constants associated with the carbon fiber oxidation and k pc , 

E pc and n pc are the Arrhenius constants associated with the pyrocarbon coating 

oxidation. The quantity p* ox is the ambient oxygen vapor pressure in the TGA 
experiments in which the values of the Arrhenius constants were determined. 

Furthermore, the quantities y/ c * and y/ pc have been introduced as the fraction of the 
fiber surface area and fraction of pyrocarbon coating surface area that is exposed to 
oxygen in the pore volume. We refer to these quantities as the exposed surface area 
fractions. The value of these parameters will also depend upon the pore morphology. 


3. FINITE ELEMENT FORMULATION 

Applying the Bubnov-Galerkin method [9], the finite element form of (2a) and 
(2b) become 


. m ox i . 

f J -dD e + f Nfr/-J ox dD e 

J * RT J dt J 

D e D e 
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(10a) 


- 1 N A^ N {*cfj + *pcj)‘^- I ^Y~~dD e =o 

D e ° D e 


and 


e & co / „ - 

f Nrf — —Ns —~J-dD e + f NjV ■ J co dD e 

J RT J dt J 1 co 

D e D e 

(10b) 

+ | N&*&LN J (* lfJ+ X pej yD'- | N$Pf«LdD°= 0 

D e ° D e 


where N t are the element shape functions and D e is the domain of each element and 
where p OXi , p coi , 5R c f. an ^ ^ pci are the element nodal values of p ox , p co , 9? c y and 
91 p C , respectively. 


Using the Product Rule of differentiation along with Gauss’ Theorem [10] and 
upon substituting equations (7a) and (7b), the second term in equations (10a) and (10b) 
can be expanded as 


\N i WJ ox dD e 

D e 


= $N,J ox -ndT e + jVN r pS x —kVp ox dD‘ 
~ ~ Pg 


D K 
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+ J VNfDj,, 

D e 


Max 

RT 


' £co\f A Vp ox dD e + f Vtf, •/>&-!_ k-Vp co dD‘ 
{P ) ^ Pi 


(11a) 


- J VNrDji 

D e 


Max 

RT 




P J 


■Vp co dD e 


and 


\N,VJ co dD e = ^N l J co ndT e + \VN r pP 0 —kVp c 0 dD e 
D e r e D e 


+ J Wi-Dm 

D e 


M, 


co 


RT 


(PoAfd.Vp dD * + f VN r pP,-Lk-V Pm dD 

{P ) ^ Pi 


- J VNi-Dm 

D e 


Mco 

RT 



\ 

<P 


P 


A 


■VPoxdD* 


(lib) 


where T e is the boundary of the element and n is the outward unit vector normal to the 
boundary. The boundary integral terms in equations (11a) and (lib) allow for the 
application of unconstrained boundary conditions. As indicated, the integration is 
performed in a closed path around the element boundary. The integral terms are only 
nonzero for elements where mass flux boundary conditions are imposed. 
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In order to linearize equations (11a) and (lib), the partial pressures which are 
operated on by the gradient operator are treated as the solution variables and it is assumed 
that the partial densities and the partial pressures in brackets in (11a) and (lib) are 
constant within the element over the duration of any arbitrary time step. Therefore, the 
partial pressures in brackets and the partial densities may be taken outside the volume 
integral and the gradient terms Vp ox and Wp co are replaced with the approximations 

and VN iPcoi . Furthermore, it is assumed that the temperature, the gas 

viscosity, the diffusivity and the volumetric porosity are all constant within the element 
(although these quantities may vary from element to element) and therefore these 
quantities may also be taken outside the volume integrals in equations (11a) and (1 lb). 

After substituting equations (11a) and (lib) into (10a) and (10b), and recognizing 
that the element nodal values of the partial pressures, p ox . and p co . , as well as the 

temporal derivatives, dp ox . / dt and dp co . / dt , are not functions of the spatial variables, 
equations (14a) and (14b) may be written in the matrix form 


COX 

0 


0 

d 

p 0Xj 

cco 

dt' 

Pco j 

* '* J 


KOX KXC 
KCX KCO 


r 

Pox j 

k 


FOX) 

FCO) 


( 12 ) 


Expressions for the elements of the capacitance and stiffness matrices and the force 
vector in equation (12) are listed in the Appendix. In the expressions in the Appendix, the 
symbol ( A ) above the quantity denotes the element average value of that quantity. 
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Equation (12) can be represented more concisely as 


[c}„^{P} n +[K]„{P} n ={F}„ 


( 13 ) 


where {p} represents the solution vector \p ox . , p COj j and where the subscript n 

indicates the matrices and vectors at time t n . Using the Backwards Difference Method 
[10], the time rate of change of the partial pressure vector is 


<*{4. Mt-W-i 

dt A t n 


where A t n = t n -t n _ j. Substituting (14), equation (13) may be rewritten in the simple 
matrix form 




(15) 


where 


k#i=TrkL+W» 


(16a) 


and 
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(16b) 


m n 

Thus, the partial pressures p ox and p co at each node are determined by the 
solution of equation (15) at each time step. The effective stiffness matrix \K e ff \ n an d 
effective force vector ^E e ff } for time t n are calculated using the expressions in the 

Appendix, equations (16a) and (16b) and the partial pressure values determined at the 
previous time step, t n _\ . A flow chart illustrating the steps involved in the numerical 
solution routine is shown in Fig. 2. 

4. APPLICATION OF THE NUMERICAL APPROACH: TGA SIMULATION 

4.1 Experiment Description 

Halbig [4] measured the weight loss behavior of 2-D plain weave C/SiC 
laminated composites. The material was fabricated by GE Power System Composites of 
Newark, DE. The material was fabricated with T300 carbon fibers and infiltrated with a 
SiC matrix through a chemical vapor infiltration (CVI) process. Specimens were 
machined into 2.54 cm-long prismatic bars with a rectangular cross-section of 0.3175 cm 
x 1.27 cm. The specimens were machined such that the through-thickness specimen 
direction (0.3175 cm) was normal to the fabric plane. The weight loss was measured in 
an environment of pure oxygen at 1 atm pressure. The residual weight was measured and 
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recorded continuously with time. A schematic of the TGA experiment and a description 
of the experimental procedure are given in Opila [11]. TGA measurements were 
performed at a number of temperatures. The results of these measurements are shown in 
Fig. 3 where the residual weight fraction (ratio of instantaneous weight to initial weight) 
is plotted versus time for each temperature. 

Fig. 3 illustrates the significant influence of temperature on the oxidation and 
weight loss behavior. Below 700 °C, the rate of oxidation and weight loss increases with 
increasing temperature. Above 700 °C, the rate of weight loss decreases with increasing 
temperature. In addition, at temperatures below 700 °C, the weight loss profile has a 
sigmoidal shape. That is, initially, the oxidation and weight loss rates are low. The rate of 
oxidation increases with time and eventually reaches a constant rate at a percent weight 
loss of approximately 15%. 

It should be noted that the sigmoidal shape of the weight loss curves at the lower 
temperatures was also observed by Lamouroux, et al. [1] while performing TGA 
experiments on T300 carbon fibers coated with a pyrocarbon coating as well as by Halbig 
[12] while measuring the weight loss behavior of T300 fibers. As Ismail [13] suggests, 
the rise in oxidation rate within the initial portion of the weight loss curve is most likely 
associated with the development of porosity in the outer surface of the carbon fibers. As 
porosity increases, more active sites become available for oxygen, thereby increasing the 
rate of carbon oxidation. The influence of porosity development on the rate of oxidation 
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is only observed in the initial portion of the weight loss curve, presumably when the outer 
portion of the fiber is being oxidized. 

At temperatures of 700 °C and above, the weight loss curves no longer exhibit the 
sigmoidal shape. The weight loss profiles maintain a constant sign of curvature (concave 
up and right). This is the case in both the C/SiC composite weight loss behavior shown in 
Fig. 3 and in the T300 weight loss behavior report by Halbig [12] and by Lamouroux, et 
al. [1]. Perhaps at the higher temperatures, the rate is no longer dependent upon the fiber 
porosity or perhaps, at the higher temperatures, the thermal expansion of the fibers results 
in sufficient porosity. 

4.2 Numerical Simulation Approach 

Two-dimensional, three-node triangular elements were chosen to implement the 
finite element method. The mesh used for the simulation of the TGA experiment, shown 
in Fig. 4, represents a cross-sectional slice of the TGA specimen. The boundary and 
initial conditions which were applied to the finite element mesh are also shown in Fig. 4. 

The material constants used for the numerical simulation are listed in Table 1. In 
this numerical simulation, we will treat the carbon fibers and pyrocarbon coating as one 

constituent. We therefore make the substitution y/ c ^ = y/ pc = ip . Furthermore, we 
assume that the Arrhenius constants are the same for both constituents. The volumetric 
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porosity was assumed to be on the order of 10% and an initial combined carbon volume 
fraction was assumed to be 50%. 

In the oxidation of carbon fibers and pyrocarbon coating in C/SiC composites, 
multiple reaction mechanisms are possible (Refs. [3], [8], [14]). For the purpose of 
demonstrating the numerical method, we will concentrate our attention on the reaction 
C + O 2 — ^ CO 2 . Thus, — &2 — 1 . 


The values of k Q and E a in Table 1 were obtained by Halbig [12] by conducting 

weight loss measurements on T300 carbon fibers. The activation energy was calculated in 
the usual manner, by plotting the natural log of the rate of weight loss versus the inverse 
of the absolute temperature. The slope of this plot yields the value of the activation 
energy. The rate of weight loss was obtained from the slope of the weight loss curves at 
the midpoint of the curves. In this manner, the variation of the rate of carbon fiber weight 
loss with temperature is input into the numerical method. 

As a result of this approach, however, the sigmoidal signature of the weight loss 
curves at the lower temperatures, that is evident in Fig. 3, will not be passed into the 
numerical solution method. Recall that the sigmoidal shape is caused by the time- 
dependent formation of porosity in the carbon fibers in the initial portion of the weight 
loss curves. Since a model which accounts for fiber porosity development and the 
associated rise in carbon oxidation rate is not part of the present numerical approach, we 
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can not expect the numerical solution to reproduce the sigmoidal shape of the weight loss 
curves at the lower temperature regime. 

Although 2-D plain weave C/SiC composites are orthotropic with regard to many 
material response properties, including the gas transport parameters of material 
permeability and areal porosity, we will assume, for the sake of simplicity, that the 
material is isotropic in regards to these transport parameters. Therefore, the permeability 
and areal porosity tensors take the form 



'o 

o 

*4e 


1 

o 

o 

1 

k = 

o 

o 

and <p A = 

o 
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o 
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1 

o 
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where k and <p A are scalar quantities representing the material permeability and areal 
porosity, respectively. 

Recognizing that the weight fraction is equivalent to the composite density 
fraction, the weight fraction remaining at each time step may be calculated with the 
numerical solution by: 1) calculating the volume average carbon volume fraction within 
the finite element mesh v c and 2) employing the rule of mixture expression to calculate 

the instantaneous composite density, p - p c v c + PsiC v SiC » where p SiC is the density of 
silicon carbide and v Si £ is the volume fraction of the matrix. We will use this approach 
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to calculate the weight fraction remaining at each time step, in order to compare the 
numerical results with the measured TGA response. 


5. NUMERICAL RESULTS 


In order to investigate the effect of material permeability on the oxidation 
behavior in the TGA experiment, a series of numerical simulations were performed using 
a wide range of permeability values. The results of these simulations revealed that the 
value of the material permeability had little or no effect on the oxidation behavior in the 
TGA experiment. The fact that material permeability has no influence on the oxidation 
behavior in the TGA experiment simulation indicates that the primary mode of oxygen 
diffusion in the TGA specimen is not Darcy flow (total gas pressure gradient-driven 
flow). This, of course, is somewhat expected, since under the conditions that have been 
imposed, the presence of significant total gas pressure gradients would not be expected. 
For all remaining numerical simulations, we will assume a value of 10" 4 m 2 /MPa-sec for 
the permeability-to-gas viscosity ratio. 

As previously discussed, the pore structure of C/SiC composites is largely 
attributed to the residual stresses in the SiC matrix that result from processing and these 
residual stresses result in matrix cracks and debonds between the carbon fibers, the 
pyrocarbon coating and the SiC matrix. This crack network in combination with the void 
volume constitutes the pore volume of the porous material. Since the magnitude of the 
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residual stresses is a function of the temperature, namely the amount of departure from 
the processing temperature, the amount of cracking and crack opening will also depend 
upon the temperature. Likewise, one would expect that the crack openings and interfacial 
debond separations will close upon reheating, and that the amount of closure is a function 
of the temperature. We also recognize that as the areal porosity and the surface area 
fractions are dependent upon the pore morphology, the value of these two material 
parameters must be affected by crack opening and closing. Specifically, larger crack 
openings and interface separations would result in a larger value for the areal porosity 
due to a smaller value for the tortuosity and, to a lesser extent, a larger value for the 
volumetric porosity. Larger crack openings and interface separations would also result in 
a larger value for the surface area fraction, since more of the carbon surfaces would be 
accessible to oxygen. Given these considerations, it is logical to assume that the areal 
porosity and the surface area fraction will be a minimum at the processing temperature 
and that the value of these parameters will increase as the temperature decreases from the 
processing temperature. 

Unfortunately, the values of these two parameters are not available for the C/SiC 
material that is the subject of this study and, in fact, it may be impossible to measure the 
in-situ value of these parameters at elevated temperatures. For this reason, a series of 
numerical experiments were performed using the numerical method and, through trial and 
error, the variation of the areal porosity and the surface area fraction with temperature 
that yielded the best match with the measured oxidation behavior was determined. This 
was based on a comparison between the rates of oxidation predicted with the numerical 
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method and the measured rates of oxidation at all TGA test temperatures. The variation 
of the areal porosity and surface area fraction with temperature that produced the best 
correlation between the measured and predicted rates of oxidation is shown in Fig. 5. We 
note that the variation of these two parameters with temperature shown in Fig. 5 is 
completely consistent with our intuition, given their dependence on the pore morphology 
and the dependence of the morphology on temperature. 

The comparison between the weight fraction remaining versus time profiles that 
were obtained with the numerical solution and the measured response is shown in Fig. 6. 
The numerical results are shown with solid lines and hollow data points and the measured 
results are shown with dashed lines and filled data points. In Fig. 6a, the results at 700 °C, 
800 °C, 900 °C and 950 °C are shown using squares, circles, triangles and diamonds, 
respectively. In Fig. 6b, the results at 600 °C and 700 °C are shown using circles and 
squares, respectively. 

Note the close agreement between the measured and predicted profiles at 800 °C, 
900 °C and 950 °C, whereas the agreement at 600 °C and 700 °C is not as close. The 
discrepancy lies in the fact that the measured weight loss curves are sigmoidal and, as 
discussed previously, the numerical method neglects to account for the initial low rate of 
carbon oxidation. We note, however, that there is close agreement between the slopes of 
the measured and predicted weight fraction curves at 600 °C and 700 °C for all points 
below a weight fraction remaining of 0.9. 
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In Fig. 7, the spatial distribution of the carbon volume fraction is plotted at 1, 2 
and 3 hrs for the numerical simulation at 700 °C and, in Fig. 8, the spatial distribution of 
the carbon volume fraction is plotted for 10, 20 and 30 hrs for the numerical simulation at 
950 °C. In addition to the fact that the rate of oxidation is much slower at 950 °C, we note 
the distinct difference between the carbon volume fraction distribution shown in Fig. 8, 
and that which is shown in Fig. 7. At 700 °C, the carbon volume fraction distribution is 
much more uniform. The carbon volume fraction profile at 950 °C, however, becomes 
more and more slender with time. Carbon oxidation occurs at the outer edge of the 
specimen with little or no oxidation in the interior. As time progresses, the delineation 
between the oxidized and non-oxidized zones marches inward from the outer surface, the 
so-called shrinking core oxidation behavior [4], It is the variation of the areal porosity 
with temperature that is responsible for the distinct difference in the oxidation patterns 
shown in Figs. 7 and 8, since areal porosity regulates the availability of oxidant in the 
interior of the specimen. 

The difference between the oxidation behavior in the two temperature regimes, 
which is illustrated in Figs. 7 and 8 has been observed by Verrilli and Calomino [15] in 
constant-load rupture test specimens tested in a partial-oxygen environment at 800 °C and 
1200 °C as well as by Halbig [16]. They observed a shrinking core oxidation pattern at 
the higher temperature and a more uniform oxidation pattern in the specimens tested at 
the lower temperature. The temperature at which the oxidation behavior transitions from 
a somewhat uniform pattern to a shrinking core pattern is somewhere between 700 °C and 
950 °C in the TGA experiment, whereas this transition temperature is between 800 °C 
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and 1200 °C in the rupture tests. The fact that the transition temperature is higher in 
rupture test specimens than in TGA test specimens makes physical sense when we 
consider that the rupture test specimens are loaded in tension and since a tension stress 
will tend to increase the areal porosity by increasing the volumetric porosity and by 
lessening the tortuosity. 


6. CONCLUDING REMARKS 


A mathematical theory and an accompanying numerical scheme have been 
developed for predicting the oxidation behavior of C/SiC composite structures. The 
theory is derived from the mechanics of the flow of ideal gases through a porous solid. 
The result of the theoretical formulation is a set of two coupled nonlinear differential 
equations written in terms of the oxidant and oxide partial pressures. The numerical 
method is based upon the solution of the two nonlinear differential equations using the 
Bubnov-Galerkin finite element method. The nonlinear differential equations are 
linearized within each time step and solved over the time domain in a piecewise linear 
manner. This is achieved by continuously updating the system stiffness matrix and 
system force vector based on the values of the solution variables determined in the 
previous time step. The end result is a numerical scheme capable of determining the 
variation of the local carbon oxidation rates as a function of space and time for any 
arbitrary C/SiC composite structure. 
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In the application section, the numerical method was applied to simulate the 
oxidation and weight loss behavior of C/SiC composite specimens in the 
thermogravimetric experiments performed by Halbig [4], To this end, it was necessary to 
determine the variation of the areal porosity and the surface area fraction with 
temperature. The numerical method was used to deduce the variation of these two 
material parameters with temperature. 

The numerical simulation method was successful in reproducing the carbon 
volume fraction spatial distribution patterns which have been observed in C/SiC stress 
rupture specimens, distribution patterns which are characteristic of the test temperature. 
Temperature was shown to have a profound influence on the oxidation behavior, both by 
its direct influence on the chemical kinetics and by its effect on the areal porosity and 
surface area fractions. 
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Appendix: Expressions for the Elements of the Capacitance and Stiffness Matrices 

and Force Vector in Equation (12) 
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Figure Captions 


Fig. 1. Optical Microscopy of 2-D C/SiC Composite (Courtesy of Southern Research 
Institute, Birmingham, AL). 

Fig. 2. Flow Chart Showing the Steps Involved in the Numerical Solution Routine. 

Fig. 3. Weight Loss Versus Time Profiles for Various Temperatures Measured in the 
TGA Experiment (Halbig [4]). 

Fig. 4. Finite Element Mesh and Boundary and Initial Conditions Used for the Numerical 
Simulation of the TGA Experiment. 

Fig. 5. Plot of the Proposed Dependence of Areal Porosity and Exposed Surface 
Area Fraction with Temperature. 

Fig. 6. Comparsion of Measured and Predicted Weight Fraction Remaining Versus Time 
Profiles for 600, 700, 900 and 950 °C. 

Fig. 7. Predicted Spatial Distribution of Carbon Volume Fraction at 1, 2 and 3 hours at 
700 °C. 

Fig. 8. Predicted Spatial Distribution of Carbon Volume Fraction at 10, 20 and 30 hours 
at 950 °C. 
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Table 1. List of Material Constant Values Used 
for TGA Oxidation Simulation 


Material Constant 

Value 

p c (g/cc) 

1.74 

=*f( sec- 1 ) 

6452.35 


1 

E* = (kJ/mole) 

118.3 

v; + v" 

c/ pc 

0.5 

<t> 

0.1 
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Fig. 2. Flow Chart Showing the Steps Involved in the Numerical Solution Routine. 
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Fig. 3. Weight Loss Versus Time Profiles for Various Temperatures Measured in the 
TGA Experiment (Halbig [4]). 
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Fig. 4. Finite Element Mesh and Boundary and Initial Conditions Used for the Numerical 
Simulation of the TGA Experiment. 


37 




Areal Porosity 
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Fig. 5. Plot of the Proposed Dependence of Areal Porosity and Exposed Surface 
Area Fraction with Temperature. 
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a. Results at 700 °C, 800 °C, 900 °C and 950 °C. 
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b. Results at 600 °C and 700 °C. 


Fig. 6. Comparsion of Measured and Predicted Weight Fraction Remaining Versus Time 
Profiles for 600, 700, 800, 900 and 950 °C. 
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